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Summary. This paper is concerned with core collapse in isolated clusters of 
identical point masses. Equations are derived which determine the evolution 
of the core which would result from given deviations from a Maxwellian 
distribution. Assuming certain forms for the deviations, the predicted 
evolution of the core is compared with the results of Monte Carlo calcula¬ 
tions. The deviations corresponding to stars moving in nearly radial orbits 
appear ultimately to dominate the evolution of the core. 

The deviations are maintained because stars of high energy cannot relax 
fast enough to thermalize with the evolving core. Under certain crude 
approximations, the deviations for stars in radial orbits can approach a 
homological form. In a better approximation, it is shown that the Fokker— 
Planck equation for the deviations admits an essentially unique homological 
solution which also determines the evolution of the core. Remaining approxi¬ 
mations prevent definitive comparison with Monte Carlo calculations, but the 
qualitative agreement is satisfactory. 

1 Introduction 

Numerical experiments apart, our understanding of the ‘gravothermal catastrophe’ is entirely 
derived from considerations of thermodynamic stability as described by Lynden-Bell & 
Wood (1968), following the work of Antonov (1962). Their conclusions were beautifully 
confirmed in numerical work by Larson (1970a), who investigated a model of self-gravitating 
TV-body systems enclosed in a reflecting sphere, as envisaged by Lynden-Bell & Wood. 

Nevertheless the relevance of the gravothermal catastrophe to models of isolated clusters 
is less clear. Though core collapse in such systems is now a well-attested fact (e.g. Aarseth, 
Henon & Wielen 1974; Spitzer & Thuan 1972), at least in the absence of energizing effects 
such as strong tidal shocks (Spitzer & Chevalier 1973) or energetic binaries (Hills 1975; 
Heggie 1975), the underlying mechanism is not clearly identified. 

In systems with stars of more than one mass, the tendency to equipartition is 
undoubtedly important (Spitzer & Shull 1975b; Aarseth 1974). In simple systems, however, 
three mechanisms have been suggested. They can be distinguished by the three different 
supposed sinks of the energy which is released by contraction of the core. One suggestion is 
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that escaping stars carry off this energy. Another is that it is absorbed in the evolution 
of a halo, i.e. the bound stars of high energy. (Throughout this paper, gravitational potential 
per unit mass is held at zero at the centre of the cluster.) Finally, in the view of Lynden-Bell 
(1968), the energy is absorbed mainly by stars which mostly occupy regions of the cluster 
where the density exceeds about 1/709 of its central value. 

The most careful discussion of numerical results with respect to these mechanisms is by 
Spitzer & Thuan (1972). They showed that, during the late stages of core collapse, the first 
two mechanisms are inadequate, and presumed that the third mechanism was likely to be 
responsible. On the other hand, Lightman & Shapiro (1978) have smudged the distinction 
between the mechanisms, arguing that the evolution of a core can be understood (and, 
indeed, predicted) on the model of escape and halo evolution. This is possible, they argue, 
provided that by ‘escape’ one understands ‘escape from the core’, and by ‘halo’ one means 
the entire region outside the contracting core. (Here the words ‘core’ and ‘halo’ have a 
spatial meaning, but for theoretical purposes a distinction in terms of energies may be more 
appropriate. Therefore in what follows we consider the core to consist of the stars with low 
energy, where the distribution is nearly Maxwellian, and the halo to consist of the stars 
of high energy, where deviations from a Maxwellian become large in some sense which we 
need not make precise.) 

Despite these differences over the understanding of the basic phenomena at work, there 
is quite close unanimity in the predicted evolution of the core. Even the crudest models 
based on escape yield the result that the central density and temperature very like specific 
powers of T = t 0 — t , where t is the time and t = t 0 at the moment of complete collapse. 
Quantitatively very similar behaviour is also predicted by a more elaborate model (von 
Hoerner 1968) based on a simplified diffusive picture of the evolution of clusters. Such 
predictions are nicely verified by detailed analysis of numerical models (Davoust 1977). 

Whatever the mechanism for core collapse, in simple systems it would not occur if the 
distribution function were precisely Maxwellian. The most rigorous contribution of the 
present paper is the derivation of a pair of equations which describe the evolution of the 
core due to deviations from a Maxwellian distribution. These equations can be used to 
discuss the influence of anisotropy on the evolution of the core, and also the effect of certain 
approximations made in Monte Carlo calculations. The evolution of the core can be 
calculated if specific forms are chosen for the deviations, and certain numerical results on 
the evolution of cores can be interpreted quite satisfactorily using these calculations. 

The two equations can be regarded as ‘equations of motion’ for the gravothermal 
catastrophe, or core collapse. On the other hand they are not closed, since the deviations 
from a Maxwellian must be specified. In the remainder of the paper, which is much less 
rigorous, we investigate how these deviations evolve. It will be shown how, on certain 
assumptions, the part of the deviations which corresponds to stars in nearly radial orbits 
tends towards a homological form, at least for high energies. Finally, using an improved but 
still very crude treatment, we show how a homological solution may lead to core evolution 
which is suggestively similar to that noted in numerical experiments. 

2 Equations for core collapse 

2.1 BASIC EQUATIONS AND ASSUMPTIONS 

Our starting point is the Boltzmann equation in which collisions are represented as in the 
Fokker-Planck equation: 

3/ 3/ 30 3/ 

3 1 dx t dx; 3 Vi 
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(Rosenbluth, McDonald & Judd 1957). Here /(x, \,t) is the number density of particles, all 
of mass m, with position vector x, and velocity v at time t; i= 1, 2, 3 labels Cartesian 
coordinates and (j) is the gravitational potential per unit mass, defined to vanish at the centre 
of the cluster and increase outward from the centre. Finally C(f, f) is the Fokker-Planck 
collision term, defined by 


C(f a ,fb )= 8trG 2 m 2 In N 


_9_/ 9 IA 1 d 2 / 
dVj \ 9u,/ 4 dVidVj \ 


9 2 I, 


9 Vj 9 Vj 


)1 


( 1 ) 


for any two distributions f a) f b , where G is the constant of gravitation, TV is the total number 
of particles, and the two ‘potentials’ are 




/»<v'wv 


v —v 



( 2 ) 


The term C(f ai f b ) describes the rate at which f a changes as a result of small-angle scattering 
off a distribution given by f b . 

If the system is large, spherically symmetric, non-rotating and approximately in a steady 
state, we can write approximately f=g(E, J , t), where E and / are, respectively, energy and 
angular momentum per unit mass. Then g satisfies the equation 


dg 3 g 30 

-+— -- =<C(/,/)> 
3 1 3 E <3 t) 


( 3 ) 


where < > denotes the time-average around an orbit of energy E and angular momentum J 
(see,e.g. Henon 1973)*. 

We assume that/is nearly Maxwellian for low values of E , and write 

g = A exp(-BE) + g 1 (E,J, 0=£o + £i (4) 

where A and B are functions of time, and gi is the deviation from a Maxwellian. It is 
specified uniquely if we require 

3gi 

£1 (0,0,0 = — (0,0, f) = 0, (5) 


which means that we choose the Maxwellian that best fits the height and slope of the correct 
distribution at zero energy. Thus, ifis neglected, the density and velocity dispersion at the 
centre of the cluster are given by 


p(0) = mA 



3/2 


and 


t4(0) = 3 B~\ 

respectively. Also, we denote the central relaxation time by 0(0). 

Let us consider qualitatively the size of g if at energies comparable with the mean kinetic 
energy per unit mass at the centre of the cluster. Depending on the starting conditions, g x 


For example, in the important case J = 0, we would have 


< 0 > = 


/' 


(pdr 


\JE - <t> 



dr 


where r is the distance from the centre of the cluster and 0(r o ) = E. 
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may be comparable with g 0 initially. However, after a few central relaxation times, colli- 
sional effects will have reduced g t greatly. Thereafter, g 0 evolves on a time-scale much greater 
than f r (0), and we shall assume that the same is true of g 1 . 

On the basis of this discussion, we shall assume that gi and its derivatives are much 
smaller than the corresponding expressions for the Maxwellian, for low values of E. Thus 
the left-hand side of (3) becomes 


) ( 6 ) 

approximately, where a dot denotes a time-derivative. Hence, provided 00/3 1) can be 
evaluated, we can find A and B (which suffice to define the evolution of the Maxwellian) 
by expanding the Fokker—Planck equation to order E , regarding E as a small quantity. In 
physical terms, we need only consider the effect of encounters on stars which remain close 
to the centre of the cluster. 

If the potential near r- 0, where r is the distance from the centre of the cluster, is 
approximately that of the isothermal model defined by g 0 , we have 


30 

A exp (-BE) -ABE exp (-BE) -AB exp (~BE)( — 

3 1 


2'nGmA/2'nf 2 

0 = —-—• l — l r 2 + higher terms. 

Hence 

1 30 _A 3 B 
0 3 1 A 2 B 

Since orbits in this potential are approximately simple harmonic, 

< 0 > = -£, 

2 


and so 



Expanding (6) to order E, we have 

A + E (-^BA - l -AB^. 

Turning now to the right-hand side of (3), we write 


(7) 


( 8 ) 


(9) 


/=/o+7i, (io) 

which is the same decomposition as in (4) except that the functions are expressed in terms 
of x, v and t. Then the collision term becomes 


< C(1 f) > = < c(/o, fo) > + < c(/o, A) > + < C(A, fo) > + < C(A, A) >, (11) 

because of its bilinear dependence on f. Here we cannot neglect the remaining terms in 
comparison with the first, since f 0 is Maxwellian and so C(/ 0 ,/o) = 0. We shall neglect the last 
term (though it can be treated in exactly the same way as the second and third terms), and 
obtain 


<C(f, /) > - < C(/ 0 , A) > + < C(f u f 0 ) >. 


( 12 ) 
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Interpreted physically, the first term on the right represents the rate at which the 
Maxwellian is destroyed by scattering off the deviations from a Maxwellian, and the second 
term represents the rate at which the deviations are being thermalized by scattering off the 
Maxwellian. 

Before we evaluate these terms, one more step is necessary. Let v t be the radial 
component of velocity at a point, and u 2 , u 3 the two tangential components. By symmetry 
about the radial direction, a Taylor expansion of /(x, v, t) in velocities will have only even 
powers of v 2 and u 3 . Hence an expansion of gi in powers of J must have only even powers of 
/, i.e. we can write 

gi(E, J, t) = g l0 (E, t) + J 2 g 10 {E, t), (13) 


if we neglect terms with higher even powers of J, where 

1 #gi 


g 10 (E, t)=gy(E, 0,0, g 20 (E , t) = 

Now since 
J 2 = r 2 (v\ + uf), 


2 dJ 


J=o 


(14) 


the neglected terms have a factor of r 4 or a higher even power of r. These factors survive all 
of the velocity integrations and differentiations in the evaluation of the collision terms, and 
so result in terms of order E 2 at least. Therefore we can evaluate (12) assuming the form 
(13) for the deviations from a Maxwellian. Finally, (5) implies that 


£io(0> t) 


dgio 


dE 


0, i.e. 


e = 0 


A 1 ^ 

gMt) -2-W 


E 2 + 0(E 2 ). 


(15) 


E= 0 


2.2 DERIVATION AND DISCUSSION OF THE EQUATIONS 

The evaluation of the collision terms for small values of | v |, corresponding to orbits of low 
energy, is a little lengthy but straightforward. The evaluation of h and / 2 is precisely 
analogous to the calculation of the electrostatic potential near the centre of a distribution 
of charge. Where a term v 2 occurs in the evaluation of C(/,/), we use the result 

(v 2 ) = E (16) 

for simple-harmonic orbits. Details of these calculations are given in the Appendix and the 
result is 

1 

r 


4 / B \ 3/2 ( r°° 

+ # gyoEdE - B~ l g 20 

E = 0 7rGmA\2iT/ \ J 0 

(18) 


i . 

-B = 

r 


1 „ r 
-B 2 
3 Jo 


- A p 10 o-l d ^0 

gl0 dE--B — 


A = -AB (* 

•V 


giodE 


(17) 
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where T= \6ir 2 G 2 m 2 In TV. For later work we note that the central relaxation time is given 
approximately by 



(19) 


where we have used equation (29) of Spitzer & Hart (1971a). 

By (1), the collision term C(f a . f b ) involves f a and integrals of f b . Therefore, in equations 
(17) and (18), terms involving integrals result from the collision term C(/ 0 , /i), and the 
remaining terms on the right-hand side arise from the term C(/i, / 0 ). Those involving g 10 
show how the isothermal core would evolve if the deviations from a Maxwellian were 
isotropic, and the other two terms show how the evolution is modified by the fact that stars 
moving on non-radial orbits may have a different distribution from those on radial orbits. 
It is not surprising that terms in g 20 do not occur in (17): the evolution of A could be 
determined by considering a particle with E = 0, i.e. at rest at the centre of the cluster, 
where the distribution is exactly isotropic. 

It may seem surprising that the evolution of the entire core can be determined by 
considering particles that remain close to the centre of the cluster. On the other hand, 
within the stated approximations, (17) and (18) are rigorous equations for the evolution 
of the best Maxwellian approximation to the correct distribution function at low energies. 
The important question really is, why is this Maxwellian also a good approximation for 
larger energies? Qualitatively this depends on the ratio of the central relaxation time to the 
time-scale for core-collapse. When this ratio is very small, as it appears to be in numerical 
experiments when core collapse is well under way, the Maxwellian approximation should be 
satisfactory for larger energies. 

In the Princeton Monte Carlo calculations (Spitzer & Hart 1971a,b; Spitzer & Thuan 
1972) it is assumed that stars scatter off a Maxwellian distribution of field stars which is 
determined essentially by the local density and velocity dispersion. Now the collision term 
(1) describes approximately the rate at which stars with distribution f a scatter off a 
population with distribution f b , and so f b is effectively isotropic in the Princeton calcula¬ 
tions. Also, f b appears in the Fokker Planck equation in the form of the integrals (2), and 
these give rise in turn to the terms containing integrals in (17) and (18). Therefore, since f b 
is assumed to be isotropic in the Princeton calculations, the term with Jg 20 EdE in (18) 
is effectively missing, and it is interesting to estimate its importance. If v 2 m and v 2 m are, 
respectively, the mean square tangential and radial velocities at a point, it can be shown from 
(4) and (13) that 


4n - 

v 2 m ~ 1 IB 



( 20 ) 


where the energy argument of g i0 and £20 is 0 + V 2 V 2 . Now the ratio of the third term on the 
right-hand side of (18) to the first term is 


6 

7r GmAB 



f 

Jo 


g 20 EdE 


f 

Jo 


tiodE 
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Table 1. Effect of terms in g i0 on core collapse. 


Model 

t/t v fr 

( 21 ) 


A2 (Spitzer & Hart 1971b) 


1.54 
+ 0.02 


4.9 

+ 0.0007 


Model 


D1 (Spitzer & Thuan 1972) 


T /*rh 

( 21 ) 


13.2 
+ 0.01 


2.8 

+ 0.0005 


0.6 

+ 0.0001 


Model 


Case 1 (Larson 1970b) 


T 

( 21 ) 


16 000 
0.03 


3.3 

0.001 


Notes: 


tjh is a reference mean relaxation time defined by Spitzer & Hart 
(1971a); r is time until core collapse and the time unit of Larson 
is 10 6 yr. 


and so, if we estimate v 2 ~ E ~ B l , we find that this is measured crudely by 

^tm ~ ^rm 1 

«&, - 1 IB Gp(0)Br 2 


( 21 ) 


where p(0) = mA(2n/B) 3/2 is the central mass-density. Expression (21), which we take as 
our measure of the relative importance of the terms in g 2 o in (18), can be understood 
roughly as the square of the ratio of the isothermal scale height to the radius where 
anisotropies become important. Estimating (21) at the point where = v 2 m , we find the 
approximate values given in Table 1. Although our estimates have no precision, it seems clear 
that the effect of g 20 is small, and becomes progressively smaller as the core collapses. Also 
Davoust’s analysis of one of Henon’s Monte Carlo models (Davoust 1977) implicitly shows 
that the radius at which anisotropy becomes substantial increases relative to the core radius. 
We conclude that the effective omission of the term with Jg 20 EdE in the Princeton Monte 
Carlo calculations does not lead to serious errors in core evolution. Furthermore, we shall 
also assume that the other term in g 20 in (18) becomes negligible. 


2.3 CORE COLLAPSE DUE TO A LOWERED MAXWELLIAN 

Let us suppose that the terms ing- 2 oin (18) are negligible. Now, although this would be true 
if the distribution were isotropic, we are not making that assumption. The function g 20 
describes the way in which the distribution of stars moving on non-radial orbits differs from 
that of stars moving radially, and what we assume is that only the latter is important as far 
as the evolution of the core is concerned. 

To proceed further, it is necessary to know something about g i 0 . Now there is some 
empirical evidence from Monte Carlo calculations (Spitzer & Thuan 1972) that g(E, 0, t) 
is approximately a ‘lowered Maxwellian’, i.e. 

g(E, 0, t) =“ A' {exp (-B'E) - exp (-£'0°°)} (E < 0«>) (22) 

together with 

g(E, 0,0 = 0 (E > 0oo), 

6 
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where A', B f are functions of time, and 

0oo = lim 0(r). 

y —► oo 


Furthermore a function g(E , /, f) of this form is an approximate solution of the Fokker— 
Planck equation (Michie 1963; King 1965; this paper, Section 3.2), and so it is at least of 
interest to compute the rate of core collapse with this assumption. 

The functions g 0 and g 10 are easily found from (4), (5), (13) and (22), which yields the 
results 


A =A' [1 - exp (~B'<f>oo)] 

B = B'[ 1 - exp (~B'(p oo)] -1 

(E < 0oo) 
(E > 0oo). 

Hence (17) and (18) yield 

—A ~ - A 2 exp (-£0oo) (1 - B<j)oo ) 

1 . 2 

— B ~ -AB exp (-Bcpoo) (9 - Bfpoo). 


and, when exp (— B’Q «,) < 1, 

U4' exp (- 5'0oo) [exp (- £'£) {1 + B'E} - 1] 
^ 10 I — A' exp (— B'E ) 


(23) 


(24) 

(25) 


Using (19), 

A 6 

t T ( 0) ex P (- B( l>~) - 0 

A yn 

B 4 

t T ( 0)—exp (— B<poo) (9 -B(p oo). 


(26) 

(27) 


Now the Princeton Monte Carlo calculations show that the central density and velocity 
dispersion tend to the approximate forms p(0) a r" 1 * 5 and i>m(0) a r -0,21 (Spitzer & Thuan 
1972), where r is the time until the instant of complete core collapse. Hence A cc T -1 - 2 
and B cc r 0,21 approximately . Thus the right- and left-hand sides of (26) and (27) can be 
calculated from data given by Spitzer and Thuan (1972), and these are compared in Fig. 1 
for their model D1. Other models from this paper give comparable results. 

In interpreting these results it should be borne in mind that they are based on the 
assumption of a lowered Maxwellian distribution. The support for this, both empirical and 
theoretical, is weak for small energies E , and therefore the value of d 2 g i 0 ldE 2 is especially 
uncertain. This only appears in the equation (18) for B but there is no reason to suppose 
that the uncertainty is insufficient to account for the fact that the observed and theoretical 
values of B differ even in sign. Even for large values of E , which contribute considerably to 
fogiodE, deviations from a lowered Maxwellian are substantial (Spitzer & Thuan 1972), 
and are in the correct sense to account for the disagreement in the values of ^ r (0)(i4/^4). 


3 A model of core collapse 

Equations (17) and (18) may be regarded as ‘equations of motion’ for core collapse. On the 
other hand they are not closed, since they do not determine g 10 and g 2 o; these are 
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- 002 - * 

14 12 10 T/ 8 6 4 2 0 

' trh 

Figure 1. Values of f r (0)04/4) and t r (0)(B/B) plotted against r/r rh for Model D1 of Spitzer & Thuan 
(1972). Dots represent values computed by assuming A <* r -1 - 2 and B^t 0 - 21 , and open circles are 
obtained from the right-hand sides of (26) and (27). 


determined by the deviations from a Maxwellian. To arrive at a better understanding of why 
cores collapse, we must ask how these deviations are maintained in the face of collisional 
effects which tend to remove them. In the present section we shall analyse this question. 


3.1 FOKKER-PLANCK EQUATION AT HIGH ENERGY 

A study of the lowered Maxwellian used in Section 2.3 shows that much of the contribution 
to the term fogiodE, which determines the rate at which A changes, comes from large 
values of E. Therefore we shall concentrate on the evolution of g for large values of E. 
Furthermore, the evidence of the Monte Carlo calculations is that g 2 o is unimportant for 
core evolution, and so we shall study (3) for / = 0 and large values of E. This means we 
consider stars travelling on large radial orbits. 

Now the collision term ) involves f 0 = A exp (-BE) multiplied by integrals of f x , 

while C(fx, f 0 ) involves /i times integrals o(f 0 . However (23) shows that/i does not decrease 
exponentially for 0 < E < <f ><*>, and so the ratio of C(/ 0 , /i) to C(f u f 0 ) becomes 
exponentially small for large values of E < 0oo. Hence (12) becomes approximately 

<C(/,/)>-<C(/ 1 ,/ 0 )>; 


and so (3) becomes 


dg dg /30\ 
bt bE\bt/ 


<C(/i,/ 0 )>. 


(28) 
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Physically, this approximation means that stars of high energy relax as though they were 
scattering off a pure Maxwellian. 

Next, stars on large radial orbits relax almost entirely when within the dense central parts 
of the cluster, and so we may approximate (2), with f b (\)=A exp [—5(0 + V 2 V 2 )], for 
very large values of v = | v |. This yields 

(2ir\ 3n 1 /2n\ 312 ( 1 ) 

exp (-50)--,, I 2 exp (— 50) | v + — j, 

neglected terms being exponentially smaller. (Incidentally, this type of approximation is 
precisely the same as that adopted by Michie in his paper of 1963.) By (1) the right-hand 
side of (28) becomes approximately 

C(/i,/o) - 8rr G 2 m 2 In N A(2nlBf 2 exp (- B<l>) (2- L \ -2— U fo (i _ 2.) 

I dVj V v ! 4 dVjdVj \ l v \ Bv f 




(29) 


before orbit-averaging. 

Analogously to (13), we can write 

fi=fio + r 2 v 2 t f 20 , (30) 


where /io(x, v, t)=g 10 (E, t) and / 20 (x, v, t) =g 2 o(E, 0- We shall now compute (C(f l0 ,f 0 )) 
and then give an argument to show that the other term in (30) may be neglected. From (29) 
we quickly find that 

C(fio, f 0 ) 47 jG 2 m 2 In N A (2ir/B) 3 ' 2 ( - —^ + I “ ex P (~ B( t >)• ( 3 0 

88 E ) v 

To find the average over one orbit, we require 



Most of the contribution to the upper integral comes from small values of r , and so we can 
set dtlv ~ drjlE approximately, and use (7). The lower integral, however, is 

{_ *— 

Jy/HE-t) 

and is contributed mostly from large values of r. Provided that most of the orbit lies within 
the region where l/i I ^ /o, we can write 


1 

0 - - In 
B 


C lnf l2 GmAr 2 
B m ~ 


(32) 


as for an isothermal sphere at large radius. (Some remarks on the validity of this approxima¬ 
tion will be found near the end of this section.) Hence 



exp (—BE/ 2) 
E 


(33) 
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Now let us return to the estimation of the contribution by the term r 2 vjf 20 in (30). 
Since we wish to evaluate C(f u f 0 ) only for radial orbits, on which v t = 0, the only 
contribution comes from the double v derivative in (1) acting on v 2 . Hence the contribution 
of the term r 2 vjf 20 to C(/i,/ 0 ) is of order 

Sir b 2 I 2 

— G 2 m 2 In Nr 2 g 20 -. (34) 

3 dvfdvj 

From the derivation of (29) and from (A 6 ), respectively, we find that 

b 2 I 2 12A (27 t/B) 3/2 [exp (- B(f>)/v] (v large) 

bvjbvj l 87 r(^ 4 / 5 ) exp (- B<p) ( v small). 


Using (7) and (32) we see that the major contribution to the orbit-average of (34) comes 
from large r, which is something of a surprise for a diffusion process in a star cluster. Using 
the form of b 2 I 2 jbVibVi for small v , therefore, we find that the orbit-average of (34) is of 
order 16/3 [Gm lnN/(27TB) 1/2 ]g 20 . Estimating ^i 0 ~ Bg i0 , g'l 0 ~ B 2 g l0 , we deduce with the 
aid of (31) and (33) that 


<C(r 2 v 2 t f M J 0 )) 

<C(/ 10 ,/o)> 


#20 

glO 



1 

-„ BE exp Qfi.BE). 

Gp (0)B 2 


(35) 


If we use (20) to estimate r 2 g 20 /Bg l0 ~(v 2 m — 2v 2 m )/(v 2 m — l/B) when the left-hand side is 
small, we find that (35) is of order (21) times BE exp (BE/ 2). If (21) is at most of the order 
10” 3 in late stages of core collapse, we deduce that it should be possible to neglect g 20 
for values of BE well below 10. In fact B(j) «>< 10 for all of the single-component Monte 
Carlo models described by Spitzer & Thuan (1972). The physical meaning of the neglect of 
terms in g 20 is that we effectively assume stars moving in radial orbits to diffuse more 
quickly in energy than in angular momentum. With this approximation, and combining 
(31) with (33), we arrive at the result 

<C(/i,/ 0 )> - 3 1/2 (2tt) 5/2 G 2 m 2 \n NAB~ 2 E~ X ex p {-BE 12) + - ^ 7 ). (36) 

l bE B bE ) 


The final step in the preliminary treatment of (28) is the evaluation of 00/30. However, 
this is comparatively straightforward, and use of (32) results in 



(37) 


neglected terms being exponentially smaller. Then, since we are dealing only with stars on 
radial orbits for which/= 0 , we can replace g in (28) by 


g*(E, t)=g(E, 0, t) =g 0 + g ro- 


(38) 


Furthermore, we can replace g l0 in (36) by g*, since 


dgo { 1 d 2 go 
bE B bE 2 
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Collecting these results, we finally obtain 


3p-* So-* 1 in [\ \ J\ 

— + — - - — BE + - = 3 :1/2 (2ir) 5 ' 2 G 2 m 2 In NAB~ 2 E~ X 
bt bE B\B\2 /A) 

+ 1 92g * ) 

+ B bE 2 ) 


exp (— BE/ 2) 



( 39 ) 


Equation (39) is an approximate equation for the evolution of the distribution of stars 
moving in large radial orbits, and it may be helpful to summarize here the assumptions 
made in its derivation. Starting from the Fokker-Planck equation, we assumed that stars 
travelling on radial orbits of large energy relax mostly when within the dense central parts 
of the cluster; thus the collision term was approximated by the form suitable for stars of 
high velocity scattering off a Maxwellian distribution. It was also effectively assumed that 
stars travelling on large radial orbits relax more rapidly in energy than in angular momentum, 
which has the effect of uncoupling the evolution of f 10 from / 20 . Finally, asymptotic 
expressions for the potential in an isothermal sphere were used for calculating 30/3r and 
orbit-averages. 

It is worth examining the last assumption here. Numerical experiments show that the 
density in the outer parts of an evolving cluster drops off more rapidly than in an isothermal 
model. However, comparison with graphs of the potential given by Spitzer &Hart (1971b) 
indicate that (32) is not greatly in error even when 0- O.80oo. Thus the use of this 
expression as an approximation for the potential seems empirically justified until energies 
close to the energy of escape are reached, and there (39) becomes invalid for another reason 
which has already been stated. 

Finally let us consider the energetics of core collapse in relation to (39). In part, the core 
evolves because its stars are involved in encounters with stars of the halo, and in principle 
it should be possible to show that the kinetic energy lost by the core in encounters is 
matched by a gain in kinetic energy by the halo. To do so it would be necessary to integrate 
(39), suitably weighted, over the region of phase space corresponding to those halo orbits 
which lead to interactions with stars of the core. Thus it is not a simple matter to investigate 
the energetics of core collapse using (39). 


3.2 A TWO-ZONE MODEL 

The right-hand side of (39) becomes negligible for large values of E. Physically, this results 
from the fact that the orbit-averaged relaxation time increases with increasing energy. For 
sufficiently large energies, say BE > y, the relaxation time is much greater than the time- 
scale for the evolution of the core*. Such stars cannot relax fast enough to remain in thermal 
equilibrium with the core. For large energies the distribution g* evolves almost entirely 
because these particles are moving through a time-dependent smooth gravitational field. 

In order to arrive at a better understanding of this process, let us consider an idealized 
model in which stars with energies less than y/B are completely relaxed, while those with 
higher energies do not relax at all. Thus 




3g* 1 d g 

+ — 


3 E B 3 E 2 


= 0 


(BE < 7) 


i(!M^}-° < " r>T) 


bg* bg * 
bt bE 

and we shall suppose that y is constant. 


(40) 


(41) 


* If an estimate of y is required, it may be obtained by equating the coefficients of dg*/bE on each side 
of (39), setting BE = 7, and using data in Fig. 1. In this way one obtains a value of about 7. 
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Not surprisingly, a lowered Maxwellian is a solution of (40), but we shall choose the 
solution 

g* = A exp (— BE) (42) 

in order to be consistent with (4) and (5). (Our Fokker—Planck equation (39) is invalid for 
low values of E , but we can look on it for the present purpose as a model equation which 
may help us to understand the full equation.) The general solution of (41), on the other 
hand, is 

g* = h[AB 1/2 exp (— BE)], (43) 

where h is an arbitrary function. The boundary conditions for h are an initial condition 
g*(E, 0 )=gg(E) (B 0 E>y) 

where go(if) is some function and B 0 is the initial value of B \and continuity with (42), i.e. 

= -4 (0 ex p (— y) 0 > 0 ). 

As our final specialisation, let us suppose that B/B 0 = (A/A 0 ) 6 , where 5 is a negative 
constant and a zero subscript denotes the initial value, and that AB 112 is an increasing 
function of time. (Both assumptions appear to be empirically justified by Monte Carlo 
calculations and those of Larson, with 5 ~ —0.2.) Then we may express the boundary 
conditions in terms of h, which determines h uniquely, and so we find that (43) becomes 


A 

$ 7 + (1 + /45) In — 
A 0 


Several features of this solution are interesting. By our previous assumptions, the quantity 
7 + (1 + V 2 b)\nA/A 0 increases with time. Therefore, at a fixed value of BE , the solution 
tends towards the upper form as t increases. Thus the influence of the initial conditions, 
which affects only the lower form, is gradually displaced towards higher values of the 
dimensionless energy BE. As far as the evolution of the core is concerned, the effect of the 
initial conditions gradually becomes negligible: straightforward calculation yields the result 

r°° A a B r°° 

g 10 dE = —exp (— 7 ) — — (1 + 345) exp (— 7 ) +— ~ I g*(E)dE 
J 0 B B B J 

y/B 0 

and the first term on the right grows faster than the others. Hence by (17) it dominates the 
growth of A . 

A second interesting feature of the solution for g* is that it tends to a Maxwellian, but 
one that is cooler than the Maxwellian for the core. Thus deviations from the core 
Maxwellian exist at high energies. Now the core is evolving as a result of the deviations from 
the core Maxwellian, and the present example shows rather clearly how these deviations may 
be maintained. At sufficiently large energies particles are effectively insulated thermally 
from the runaway core, at least on the time-scale of its evolution. 


BE + Vidy] 
A exp {-' 

' V28 J 


r* = 


1 + 


1 + US A 

-In 

Bq A , 


for BE 


B \ 

0 B 0 / 
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A final remark shows us the way to further progress. The solution to which g* tends in 
the above crude model is homological, i.e. it is of the form 


g*(E,t)=Ay(BE), 


where y is some function. 


(44) 


3.3 HOMOLOGICAL CORE COLLAPSE 

In a rather trivial sense, the evolution of the core must be homological if the core is approxi¬ 
mately isothermal, for different isothermal models are related by simple scalings. However 
(44) is a stronger assumption, for it means that also the deviations from the core Maxwellian 
scale in the same way. This is the sense in which we use the word ‘homological’ here. As we 
shall see, a consequence of this assumption is that parameters for the core evolve as simple 
powers of r, the time until complete collapse. 

The conjecture that core collapse is approximately homological (Lynden-Bell 1975) 
seems well supported by some aspects of numerical work on clusters. The approximate 
power law dependence on time of the central density and velocity dispersion (Larson 
1970a, b; Spitzer & Thuan 1972; Davoust 1977) is evidence for this idea, and also the 
evolution of the density profile can look temptingly self-similar (Larson 1970a). On the 
other hand it is clear from Table 1 and Davoust’s analysis (Davoust 1977) that the profile 
of the velocity dispersions does not partake of this evolution. 

We shall now search for self-similar solutions of (39), which means only that we suppose 
the evolution of g* to be self similar, i.e. the evolution of the distribution of stars moving 
in radial orbits. Thus we avoid the problem of imposing self similarity on the anisotropy 
also. 

By substituting (44) into (39) we have 


y +y'; j- — +— J = 3 y2 (2n) 5l2 G 2 m 2 In NAx 1 exp (- V 2 x)(y" +y' 


(45) 


where x =BE . A homologous solution will be possible if A/ A 2 and B/BA are constants, and 
so we suppose that 


— =3 V2 (27i) 5/2 G 2 m 2 \nNa 
A 2 


(46) 


and 

B A 
-=8 — 

B A 

where a and 6 are constants. Thus (45) yields 
y" + y f = ax exp 0 6.x){y +y'( 1 + V 2 S)}. 

By analogy with (38) we have 
y = yo + yio = exp (-x)+y 10, 


(47) 


(48) 


( 49 ) 
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where y 10 is the scaled deviation from a Maxwellian for stars on radial orbits, i.e. g xo = 
Ay X0 (BE). Hence, combining (17), (18), (46) and (47), we have 


/* oo 

2x/2/37r y l0 dx 

(50) 

Jo 


2 ( 16 yl o(0) 

3 3 r°° 

(51) 

y 10 dx 



o 


provided that we ignore terms involving g 2 o Note that the assumption of self similarity 
has proved self consistent: to obtain self-similar solutions we had to assume that a and 6 
were constants; now we see that they really are constants. 

We still do not quite have a well-defined problem, since (48) should only be valid for 
values of x well above 1. However, for values of x < 1, we know from (15) and (49) that 


y 


exp (~x) + -yio(0)x 2 . 

2 


(52) 


Therefore, in the absence of the correct equation fory for all x, the best we can do is to 
match the solution of (48) on to (52) at a value x = x m , say, where x m is of order unity. 

The problem thus reduces to the following. Taking trial values for a and 6, we integrate 
(48) from some large value of x, choosing the initial condition from the asymptotic form of 
its solution, i.e. y ~ C exp — [x/(l + Vid)], where C is a constant. At x = x m we choose C 
and yi'o(O) so as to make y and y match with (52). Finally, we compute a and 6 from (50) 
and (51), and iterate until the guessed and computed values agree. The result of this 
calculation depends on the choice ofx m , as we see in Fig. 2, which we shall discuss shortly. 

Apart from the values of a and 5, the relevance of the homological model to Monte Carlo 
a nd hy drodynamic calculations can be tested by the assumed constancy of a. Now a = 
\j2j21t x (0)(4/4), and so we already have some information on this in Fig. 1. (For this 
purpose the dots provide the relevant information.) Furthermore, if a is constant, (46) 
implies that A a 7^, where ft = — 1 exactly. Values of j3 derived from a selection of numerical 


Table 2. Values of 0. 



Model 

& 

Larson (1970a) 

Enclosed system 

-1.00 

Larson (1970b) 

Case 1 

-1.08 


Case 4 

-1.09 


Galactic cluster 

-1.24 


Galactic nucleus 

-1.02 

Spitzer & Thuan (1972) 

D 

-1.24 


E 

-1.49 


F 

-1.15 

Henon 

72 

—1.55 

Spitzer & Shull (1975a) 

K 

-1.61 


Notes: 

The value for the last five models has been obtained from the 
time-dependence of the central velocity dispersion and density 
as tabulated by Davoust (1977). For the last two models 
Davoust gives two determinations each, and these were averaged 
to yield the values listed here. 
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Figure 2. Values of a and 6 for the homological model described in the text and for various numerical 
models of clusters. Those determined from the homological model depend on x m , the value of x at which 
matching is performed; these values are plotted along the curve. The values for the numerical models 
depend on time (see Fig. 1), and the values plotted are those at the last time for which details are 
published. They are identified as follows. ES: enclosed system (Larson 1970a); I: globular cluster, case 1; 
IV: globular cluster, case 4; GC: galactic cluster; GN: galactic nucleus (Larson 1970b). Dl, El, F2 are 
models of Spitzer & Thuan (1972). 


results are given in Table 2. From these results it already seems that homological solutions 
may be of limited applicability, and several possible reasons could be suggested for this. 
Quite apart from the possibility that homological behaviour may not be inevitable, we may 
still be witnessing the effects of initial conditions. Possibly our approximation (32) for the 
potential may be insufficiently good for particles of large energy whose relaxation times 
much exceed the time-scale for evolution of the core. However, the fact that j3 =£ — 1 in most 
cases has nothing to do with the crude matching procedure adopted; the r _1 -dependence of 
A is a general prediction for any homological solution. 

Now let us consider the values of a and 6, plotted in Fig. 2. Because of the dependence 
on x m , one can make no definitive observations. On the other hand it is no discouragement 
to note that there is a reasonable choice ofx m for which both 6 and a fall within the ranges 
obtained in numerical work. 

Further differences between the homological solution derived here and the Princeton 
Monte Carlo calculations are shown in Fig. 3. It illustrates that, in the homological solution, 
deviations from a Maxwellian become important relative to the core Maxwellian only for 
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Figure 3. The function -y i0 (x) for the solution of (48) and (52) corresponding tox m = 2. Also shown 
is the Maxwellian y Q (x) = exp (-x), and the function -y 10 (x) which is the scaled version of the 
deviation corresponding to a lowered Maxwellian with B<f >« = 9.5; this is calculated from (23) and is 
indicated by a dashed curve. 


large values of x=BE , i.e. x > 12. On the other hand even at the late stages of core collapse 
in the models of Spitzer & Thuan (1972), the escape energy corresponds to x ~ 9.5, and 
deviations from a Maxwellian must be important even at scaled energies slightly lower than 
this. This is illustrated by the graph of deviations corresponding to a lowered Maxwellian. 
The remarks made in connection with the time-dependence of A are equally relevant to the 
results of the last two paragraphs. 

4 Discussion 

4.1 THE EVAPORATIVE MODEL 

Until now, the most successful quantitative treatment of core collapse has been the one 
described in its most updated version by Lightman & Shapiro (1978). It cuts through the 
complications of the Fokker—Planck equation by ignoring this and reverting to a more 
primitive model for cluster evolution. In this model, stars escape at a certain fractional rate 
per relaxation time, carrying off a small, perhaps negligible, flux of energy, and so the time- 
dependence of the mass and total energy may be found. Originally applied to escape from a 
cluster (Spitzer 1940), it is now applied to escape from the core of a cluster into the halo. 
Furthermore, in the spatial sense of the word, the ‘halo’ is now to be understood as a region 
which extends closer and closer to the centre of the cluster as the core radius shrinks 
(Lightman & Shapiro 1978). 

The mathematical expression of the evaporative model is summarized by the equations 

l^ c = _i 

N c dt T 
1 dE c p 
E c dt T 
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where /V,_ /-; c are, respectively, the number of stars in the core and its energy, T is the central 
relaxation time, and p, q are constants (von Hoerner 1976). If we estimate the parameters 
for the core Maxwellian from 


A « N c ~ 1/2 Ei n , B oc N c E~ l 

(which follow using expressions previously given for the central density and velocity 
dispersion), we deduce that 


1 dA Iq + 3 p 
A dt~ 2 T 

1 dB p+q 

B dt T 


(53) 

(54) 


Using the model of core collapse developed in Section 3.3, very similar results may be 
found. In fact from (19), (46) and (47) we obtain 


1 dA_ 1 127 

A dt' t T (0 ) V 2 “ 


1 dB 
B dt 


1 

MP) 



which are of essentially the same form as (53) and (54). However, essential differences 
appear in the way in which the constants on the right-hand sides are obtained. In the model 
discussed in this paper, they can be written in terms of the deviations from a Maxwellian 
(equations (50) and (51)), and are determined along with those deviations by solution of a 
differential equation. On the evaporative model, the determination of the constants is some¬ 
what problematic. 

On the evaporative model, p and q are sometimes calculated by assuming that, in each 
relaxation time, all stars in the tail of a Maxwellian distribution above some velocity will 
escape. Difficulties with this method include the choice of the escape velocity and 
justification for the use of a Maxwellian (von Hoerner 1976). As a modification, it may be 
noted that stars escape with very little excess energy in very large systems, and so one may 
take p = 0 (Lightman & Shapiro 1978). However, this ignores the fact that the energy of 
the core also changes as a result of encounters with stars of the halo, even when such 
encounters do not lead to escape. Finally, it is possible to determine p and q by comparison 
with numerical experiments such as Monte Carlo calculations. However, one cannot then 
have confidence that the same values are applicable under different circumstances, for 
example in the presence of large numbers of binaries. 


5 Conclusions 

In the present treatment, core collapse is regarded as the balance between two processes. 
Firstly, the core will evolve in the presence of deviations from a Maxwellian distribution, 
in accordance with equations (17) and (18). Secondly, the deviations evolve in response to 
the collapse of the core. Stars with small binding energy can relax sufficiently rapidly to 
keep the deviations small, but those with large energy cannot and so the deviations are 
maintained. Collapse of the core occurs in such a way that these two processes occur at a 
self-consistent rate. 
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A deeper reason for the fact that such a situation is possible may, perhaps, be sought in a 
discussion of the gravothermal catastrophe. If the deviations from a Maxwellian become 
important at too low an energy, they may simply thermalize away, just as an isothermal 
system enclosed in a sphere of small radius shows no tendency to a gravothermal runaway 
(Lynden-Bell & Wood 1968). On the other hand, if the deviations become important only at 
a very great energy, where the relaxation time is extremely long, core collapse will occur on 
an extremely long time-scale. It seems reasonable to suppose that, in a typical system, the 
deviations will tend to become important at the lowest energy consistent with a gravo¬ 
thermal runaway, since core evolution then takes place on the fastest possible time-scale. 

The present treatment shows that there are situations in which the initial conditions 
gradually become less important as the core evolves, and this can lead to self-similar 
behaviour. Unfortunately, most of the analysis in this paper depends on the assumption, 
supported only by numerical results, that the dependence of the distribution on/for small 
values of J is negligible in certain respects. Furthermore, only a better treatment of the 
Fokker-Planck equation for the deviations from a Maxwellian will allow us to decide 
whether the crude matching procedure we have adopted accounts for many of the 
discrepancies with the results of Monte Carlo calculations. Finally, the present treatment 
deals only with the case of equal masses, and new features might arise in more realistic 
systems (cf lightman & Fall 1978). Nevertheless, it may represent the first step in a new 
quantitative description of core collapse based firmly on the Fokker-Planck equation. 
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Appendix: evaluation of collision terms for small energies 

In this appendix we shall evaluate the right-hand side of (12) for small energies,#, as far as 
terms which are linear in E. As explained in the text, higher terms omitted in the expansion 
(13), i.e. (30), can be neglected, and so 

<C(/ /)> “ <C(/o,/io)> + <r 2 C(fo, v 2 tf 20 )) + <C(/ 10 ,/o)> + <r 2 C(v 2 t f 2(h f 0 )). (Al) 

A little preliminary reduction is useful, however. Since V$I 2 = 2 I x by (2), then (1) becomes 


C{f a J b ) = lit G 2 m 2 In N (- f a V 4 v h + ~~ ~~ 

1 uVfuVj 9 Vf 9 Vj 

Now we deal with successive terms of (Al) in turn. 


(A2) 


1 <C(/o,/ 10 )>. 

We must evaluate I 2 to terms in u 6 , since then the terms in (A2) will be correct to order v 2 , 
i.e. 0(E). Since f 10 is isotropic, we have 


(A3) 


47 T r v V 47T 

h = — /lo(v') — (v 2 + 3 v 2 )du + — /io(v')i/(u 2 + 3 v 2 )dv 

3 Jo ” 3 J y 

exactly, where some of the arguments of f i0 have been suppressed. Now/i 0 = 0(# 2 ) for small 
energies, by (15); hence the first term is 0(# 4 ), i.e. 0(u 8 ), and so it can be neglected. 
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Core collapse in clusters 175 

Furthermore the second integral can be extended to 0 < v < 00 with comparable error. Thus 

d% Sn r°° 

fio(y')v'dv' 

oVidVj 3 J 0 

and so the curly bracket in (A2), with f a replaced by f 0 =A exp {- QAv 2 + (j>)B }, comes to 
- 87 tBA f g 10 (E')dE' (1 -B<p - V 2 Bv 2 - VsBv 2 ) 


where we have written /io(V) =gio(E ) and have neglected terms of order E 2 or higher. 
For the same reason, the integral can be extended to 0 < E’ < °°. Performing the orbit- 
average using ( 8 ) and (16), we find 


<C(/ 0 ,/io)>- - TABil-ysBE) f gi 0 (E')dE f 

J 0 

if we neglect terms of 0(E 2 ), where T = 167 r 2 G 2 m 2 In N. 


(A4) 


2 <C(/ 0 ,r 2 u?/ 20 )>. 


The presence of the factor r 2 , which is of order E, means that we need only evaluate / 2 to 
order v 4 to retain all terms in (A2) of order E. Now 

I 2 = r 2 ( |v-v'|u; 2 / 2 o(vVV 

J 

in this case, and, as before, we may separate the integral into the range v^v. The 
contribution from v < v is negligible to the order in which we are interested. For v > v 
we may expand | v — v' | as a power series in v as far as terms of order v 4 , omitting terms 
which are odd in V since the rest of the integrand is even in v'. (The neglected higher terms 
would give contributions at most of order v 6 In v.) Thus 

, ( 1 y 2 1 (v • v ') 2 1 v 4 3 (v • v') 2 u 2 

l 2 v 2 2 v 4 8 v 4 4 i / 6 

5 ( V ' V ') 4 | , 2f ( > \ j3 - 

~8~^-) Vtf20(y)d 

Again, however, we can extend the integration to 0<i/<°° with an error which does not 
concern us. Now integration in spherical polar coordinates yields 



h ~ 


8tt r 


r > 5 , f,, 4v r + 3vf, - 2 v 4 -vv +v ? 

J 0 1 10 u 2 70u' 4 


dv 


where t> r , v t are the radial and transverse components of v. 

To proceed further it is convenient to choose a specific local Cartesian coordinate system 
as in Section 2.1, and so we find that the curly bracket in (A2) is 


167T 

3 


r 2 AB 


P : 

2 0 


'fiodv. 


plus terms of order E 2 . Next, writing /mIV) =g 2 o(E') where E' = <f> + Viv' 2 , we can write 
this as 

32tt 

3 


Ab( (E'-<t>)g M dE', 

J (h 
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176 D.C.Heggie 

and here we can take 0 = 0 if we want to retain only those terms of order E. Combining with 
the factor in (A2), and finding the orbit-average of r 2 by combining (7) and (8), we finally 
obtain the result 


<C(/o,^^/2o)>- -—BE 
nGm 


01 ? 


giodE', 


(A5) 


COMPLETION OF THE CALCULATION 

The remaining terms of (Al) are much simpler to compute. Since 

3 2 £io 

E= 0 


fio- V2E : 


dE 


for small E, we need I 2 only as far as terms of order v 2 . Using (A3), with/ 10 replaced by 
fo=A exp (—BE), we find that the first term is negligible again, and that the range of 
integration in the second can be extended to 0 < v < °°. Thus 


a 2 /. 


87r A 

- exp (- B(j))6 if 


3 VidVj 3 B 

whence 

4T A a 2 g 10 


(A6) 


and 


e = o 


T /B\^ 2 E 


<C(r 2 v 2 f 20 ,f 0 ))-—- - — ^2o(0). 

nGmUn' B 

Combining (A4) and (A5) with these two results, and equating with (9), we have 


A +E 


/ 3 . 

1 A 1 / 

4 \ r 

\--BA 

--AB\ =r \-AB 1 

—BE) 

\ 2 

4 / \ 

3 JJ( 


gio(E')dE r 


BE IB 


ttGm 


(Dl? 


g ^E)„E*-- — 


1 1‘ 

B\ 3/2 E ) 

+ —-(- 

- ) -^ 20 ( 0 ) 

o nGm \2 

.nf B 1 


Equating like powers of#, and simplifying slightly, we arrive at (17) and (18). 

It is to be noted that equations (17) and (18) are as rigorous as the Fokker—Planck 
equation, except for two inessential approximations: the approximation of the central 
density used to arrive at (7), and the neglect of the last term in (11). Furthermore, the 
calculation can be extended very easily to take account of the presence of different stellar 
masses. 
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